function DecodingProbability_v7()

% clear workspace
clear all
close all
clc

T=30; % transmitted packets
q=2; % field cardinality
g=[0.7 0.3]; %  layer selection probability
k1=10;
k2=10;
K1=k1;
K2=k1+k2;


% % Decoding 1. Layer(working?)
% l1_prob=zeros(T,1);
% for pkts_recv = 1:T % for each nb of transmitted packets
%     bino_vector=binopdf(1:pkts_recv,pkts_recv,g(1)); % binomial vector with the probability of 1,2,3,...,pkts_recv layer 1
%     for outcome=1:length(bino_vector) % for each outcome of the received packets
%         if outcome>=k1
%             l1_prob(pkts_recv)=l1_prob(pkts_recv)+bino_vector(outcome)*ProbMatricesWithRank(outcome,k1,k1,q);
%         end
%     end
% end



% Decoding 1. Layer(working?)
l1_prob=zeros(T,1);
for pkts_recv = 1:T % for each nb of transmitted packets
    l1_prob(pkts_recv)=MatrixFun(pkts_recv,k1,k1,g(1),q);
end

l1_prob

% Decoding 2. layer
l2_prob=zeros(T,1);
for pkts_recv = 1:T % for each nb of transmitted packets
    
    for i=0:K1
        
        if pkts_recv<K2
            l2_prob(pkts_recv)=0;
        else
            
            %             if pkts_recv<i
            %                 val1=0;
            %             else
            %                 val1=MatrixFun(pkts_recv,k1,i,g(1),q);
            %             end
            %
            %             if pkts_recv<K2-i
            %                 val2=0;
            %             else
            %                 val2=MatrixFun(pkts_recv,K2-i,K2-i,g(2),q);
            %             end
            %
            %             l2_prob(pkts_recv)=l2_prob(pkts_recv)+val1*val2;
            val1=MatrixFun(pkts_recv,k1,i,g(1),q);
            val2=MatrixFun(pkts_recv,K2-i,K2-i,g(2),q);
            
            l2_prob(pkts_recv)=l2_prob(pkts_recv)+val1*val2;
            
            
        end
    end
end

l2_prob



% Plotting
figure(1)
hold('on')
% plot(1:T,p_decode_l1)
plot(1:T,l1_prob)
plot(1:T,l2_prob)
% plot(1:T,p_decode_l3)
hold('off')
grid('on')
pbaspect([2.5 1 1])
set(gca,'XTick',0:10:T)
xlim([0 T])
ylim([0 1])

% % Save plot
% print(gcf,'uep_ew_analytic.eps')

end


function Pr = MatrixFun(pkts_recv,layer_length,rank,g,q)

Pr=0;

if pkts_recv>=rank % if we have less pkt's than required rank, no need to calculate -> impossible
    
    bino_vector=binopdf(1:pkts_recv,pkts_recv,g); % binomial vector with the probability of 1,2,3,...,pkts_recv layer 1
    for outcome=1:length(bino_vector) % for each outcome of the received packets
        if outcome>=rank
            Pr=Pr+bino_vector(outcome)*ProbMatricesWithRank(outcome,layer_length,rank,q);
        end
    end
    
else
    Pr=0;
end
% Pr;


end

% Should work
function PMWR = ProbMatricesWithRank(m,n,r,q)

% Get first set of gaussian coefficients
gc=gausscoeffs(n,r,q);

% Calculate "sum"
val=0;
for k=0:r
    % This should be the one!
    val=val+((-1)^(r-k)*gausscoeffs(r,k,q)*q^(m*k+binomcoeffs(r-k,2)-n*m));
end

% Return probability of matrix 'm'x'n' with rank 'r'
PMWR=gc*val;

end

% Should work (Tested! see bottom)
function GC = gausscoeffs(m,r,q)
if r==0
    % disp('r = 0 in gauss coeffs')
    GC=1;
elseif r>0
    % disp('r > 0 in gauss coeffs')
    
    % Calculate numerator
    num=1;
    for w=m:-1:m-r+1
        num=num*(q^w-1);
    end
    
    % Calculate denominator
    denom=1;
    for w=r:-1:1
        denom=denom*(q^w-1);
    end
    
    % Calculate gaussian coefficient
    GC=num/denom;
    
elseif r<0
    disp('r < 0 error in gausscoeffs!!!')
end


end

% Not tested!
function bc = binomcoeffs(a,k)
% As on page 123 in "A course in combinatorics"

tmp_vector=ones(2,1);
tmp_index=1;

for w=0:-1:-k+1
    tmp_vector(tmp_index)=(a+w);
    tmp_index=tmp_index+1;
end

num=prod(tmp_vector);

denom=factorial(k);
bc=num/denom;

end



